Setup

knitr::opts_chunk$set(cache = TRUE)
library(tidyverse)
library(leaflet)
library(sf)
library(dbscan)
library(spatstat)
library(splancs)
library(htmlwidgets)
library(webshot)
data_dir <- "../clean_data/"
maps_dir <- "../maps/"

EDA & Data Wrangling

starbucks <- read_csv(paste0(data_dir, "starbucks_shanghai.csv")) |>
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326) |>
  select(-city)
## Rows: 754 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): city
## dbl (2): latitude, longitude
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
starbucks
## Simple feature collection with 754 features and 0 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 121.0167 ymin: 30.73263 xmax: 121.9256 ymax: 31.62088
## Geodetic CRS:  WGS 84
## # A tibble: 754 × 1
##               geometry
##            <POINT [°]>
##  1 (121.3965 31.62088)
##  2 (121.2479 31.39764)
##  3 (121.2374 31.39053)
##  4 (121.3506 31.40225)
##  5 (121.2456 31.38443)
##  6 (121.2562 31.38369)
##  7 (121.2788 31.38513)
##  8 (121.3364 31.39046)
##  9 (121.3583 31.38815)
## 10 (121.2587 31.35715)
## # ℹ 744 more rows
shanghai <- st_read(paste0(maps_dir, "2019年县级.shp")) |>
  filter(NAME_1 == "Shanghai") |>
  select(ENG_NAME) |>
  rename(District = ENG_NAME) |>
  mutate(District = case_when(
    District == "Pudongxin" ~ "Pudong New Area",
    District == "Jingan" ~ "Jing'an",
    TRUE ~ District
  ))
## Reading layer `2019年县级' from data source 
##   `/Users/bulldf/University of Toronto/2024 Fall/STA465/project/maps/2019年县级.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2856 features and 29 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 73.50114 ymin: 3.837902 xmax: 135.0885 ymax: 53.5609
## Geodetic CRS:  WGS 84
shanghai
## Simple feature collection with 16 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 120.8523 ymin: 30.6917 xmax: 122.1129 ymax: 31.87255
## Geodetic CRS:  WGS 84
## First 10 features:
##           District                       geometry
## 1        Chongming MULTIPOLYGON (((121.4272 31...
## 2         Fengxian MULTIPOLYGON (((121.5661 31...
## 3          Hongkou MULTIPOLYGON (((121.4808 31...
## 4          Huangpu MULTIPOLYGON (((121.4904 31...
## 5          Jiading MULTIPOLYGON (((121.3318 31...
## 6          Jinshan MULTIPOLYGON (((121.0238 30...
## 7          Jing'an MULTIPOLYGON (((121.4638 31...
## 8          Minhang MULTIPOLYGON (((121.3336 31...
## 9  Pudong New Area MULTIPOLYGON (((121.8514 31...
## 10          Qingpu MULTIPOLYGON (((121.1492 31...
demographics <- read_csv(paste0(data_dir, "shanghai_2019_en.csv"))
## Rows: 16 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): District
## dbl (3): Area, Population, Density
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
demographics
## # A tibble: 16 × 4
##    District          Area Population Density
##    <chr>            <dbl>      <dbl>   <dbl>
##  1 Pudong New Area 1210.       557.     4599
##  2 Huangpu           20.5       65.1   31808
##  3 Xuhui             54.8      109.    19989
##  4 Changning         38.3       69.4   18110
##  5 Jing'an           36.9      106.    28680
##  6 Putuo             54.8      128.    23268
##  7 Hongkou           23.5       79.4   33816
##  8 Yangpu            60.7      130.    21487
##  9 Minhang          371.       255.     6876
## 10 Baoshan          271.       204.     7544
## 11 Jiading          464.       160.     3438
## 12 Jinshan          586.        80.7    1377
## 13 Songjiang        606.       177.     2926
## 14 Qingpu           670.       123.     1840
## 15 Fengxian         687.       116.     1684
## 16 Chongming       1185.        68.4     577
leaflet() |>
  addProviderTiles("CartoDB.Positron") |>
  addCircles(
    data = starbucks,
    color = "darkgreen"
  ) |>
  addPolygons(
    data = shanghai,
    color = "#4B4B4B",
    weight = 1.5,
    opacity = 1,
    label = ~District
  )
wgs_coords <- st_coordinates(starbucks)

utm_coords <- starbucks |> st_transform(crs = 32651) |> st_coordinates()

df <- starbucks |>
  mutate(
    lon = wgs_coords[, 1],
    lat = wgs_coords[, 2],
    x = utm_coords[, 1],
    y = utm_coords[, 2]
  ) |>
  st_drop_geometry()
df
## # A tibble: 754 × 4
##      lon   lat       x        y
##  * <dbl> <dbl>   <dbl>    <dbl>
##  1  121.  31.6 347914. 3499530.
##  2  121.  31.4 333421. 3474999.
##  3  121.  31.4 332409. 3474227.
##  4  121.  31.4 343190. 3475359.
##  5  121.  31.4 333176. 3473538.
##  6  121.  31.4 334182. 3473440.
##  7  121.  31.4 336336. 3473566.
##  8  121.  31.4 341819. 3474072.
##  9  121.  31.4 343900. 3473785.
## 10  121.  31.4 334376. 3470495.
## # ℹ 744 more rows
window <- convexhull.xy(df$x, df$y)
window
## window: polygonal boundary
## enclosing rectangle: [310427.1, 397312.9] x [3401137, 3499530] units
starbucks_ppp <- ppp(df$x, df$y, window = window)
starbucks_ppp
## Planar point pattern: 754 points
## window: polygonal boundary
## enclosing rectangle: [310427.1, 397312.9] x [3401137, 3499530] units
starbucks_pts <- as.points(df$x, df$y)

plot(starbucks_ppp, main = "Point Pattern of Starbucks Stores in Shanghai")

Homogeneous PPP

hppp <- rpoispp(lambda = nrow(df) / area(window), win = window)
hppp
## Planar point pattern: 736 points
## window: polygonal boundary
## enclosing rectangle: [310427.1, 397312.9] x [3401137, 3499530] units
plot(hppp)

Testing for CSR

Ripley’s K

min_x <- min(df$x)
max_x <- max(df$x)
min_y <- min(df$y)
max_y <- max(df$y)
poly <- as.points(c(min_x, max_x, max_x, min_x), c(min_y, min_y, max_y, max_y))

starbucks_seq <- seq(0, 80000, 2000)

khat <- khat(starbucks_pts, poly, starbucks_seq)
khat
##  [1]          0  242653037  790222405 1510179002 2284491196 3072215026
##  [7] 3797145392 4442974302 5023099746 5524568337 5952076682 6312154056
## [13] 6623368200 6893070622 7142798921 7378406327 7608843623 7840342986
## [19] 8048979070 8217006776 8363620828 8489521559 8598801337 8688037900
## [25] 8780845469 8875149534 8967234186 9052615422 9174949555 9309131136
## [31] 9409884731 9489957803 9557019064 9609670451 9651253736 9678783948
## [37] 9694863508 9709981396 9718607706 9723013929 9729432091
ul_khat <- Kenv.csr(length(df$x), poly, nsim = 99, starbucks_seq)
## Doing simulation  1 
## Doing simulation  2 
## Doing simulation  3 
## Doing simulation  4 
## Doing simulation  5 
## Doing simulation  6 
## Doing simulation  7 
## Doing simulation  8 
## Doing simulation  9 
## Doing simulation  10 
## Doing simulation  11 
## Doing simulation  12 
## Doing simulation  13 
## Doing simulation  14 
## Doing simulation  15 
## Doing simulation  16 
## Doing simulation  17 
## Doing simulation  18 
## Doing simulation  19 
## Doing simulation  20 
## Doing simulation  21 
## Doing simulation  22 
## Doing simulation  23 
## Doing simulation  24 
## Doing simulation  25 
## Doing simulation  26 
## Doing simulation  27 
## Doing simulation  28 
## Doing simulation  29 
## Doing simulation  30 
## Doing simulation  31 
## Doing simulation  32 
## Doing simulation  33 
## Doing simulation  34 
## Doing simulation  35 
## Doing simulation  36 
## Doing simulation  37 
## Doing simulation  38 
## Doing simulation  39 
## Doing simulation  40 
## Doing simulation  41 
## Doing simulation  42 
## Doing simulation  43 
## Doing simulation  44 
## Doing simulation  45 
## Doing simulation  46 
## Doing simulation  47 
## Doing simulation  48 
## Doing simulation  49 
## Doing simulation  50 
## Doing simulation  51 
## Doing simulation  52 
## Doing simulation  53 
## Doing simulation  54 
## Doing simulation  55 
## Doing simulation  56 
## Doing simulation  57 
## Doing simulation  58 
## Doing simulation  59 
## Doing simulation  60 
## Doing simulation  61 
## Doing simulation  62 
## Doing simulation  63 
## Doing simulation  64 
## Doing simulation  65 
## Doing simulation  66 
## Doing simulation  67 
## Doing simulation  68 
## Doing simulation  69 
## Doing simulation  70 
## Doing simulation  71 
## Doing simulation  72 
## Doing simulation  73 
## Doing simulation  74 
## Doing simulation  75 
## Doing simulation  76 
## Doing simulation  77 
## Doing simulation  78 
## Doing simulation  79 
## Doing simulation  80 
## Doing simulation  81 
## Doing simulation  82 
## Doing simulation  83 
## Doing simulation  84 
## Doing simulation  85 
## Doing simulation  86 
## Doing simulation  87 
## Doing simulation  88 
## Doing simulation  89 
## Doing simulation  90 
## Doing simulation  91 
## Doing simulation  92 
## Doing simulation  93 
## Doing simulation  94 
## Doing simulation  95 
## Doing simulation  96 
## Doing simulation  97 
## Doing simulation  98 
## Doing simulation  99
plot(
  starbucks_seq,
  khat - pi * starbucks_seq^2,
  type = "l",
  xlab = "Distance",
  ylab = "Estimated K - pi * h^2",
  main = "Ripley's K"
)

# plot upper bound
lines(starbucks_seq, ul_khat$upper - pi * starbucks_seq^2, lty = 2)

# plot lower bound
lines(starbucks_seq, ul_khat$lower - pi * starbucks_seq^2, lty = 2)

l <- function(k, h) {
  sqrt(k / pi) - h
}

plot(
  starbucks_seq,
  l(khat, starbucks_seq),
  type = "l",
  xlab = "Distance",
  ylab = "Estimated L"
)

# plot upper bound of Lhat
lines(starbucks_seq, l(ul_khat$upper, starbucks_seq), lty = 2)

# plot lower bound of Lhat
lines(starbucks_seq, l(ul_khat$lower, starbucks_seq), lty = 2)

K-S Test

ks_x <- cdf.test(starbucks_ppp, test = "ks", "x")
plot(ks_x)

ks_y <- cdf.test(starbucks_ppp, test = "ks", "y")
## Warning in ks.test.default(U, "punif", ...): ties should not be present for the
## Kolmogorov-Smirnov test
plot(ks_y)

G-function

plot(envelope(starbucks_ppp, Gest), main = "G-Function")
## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
## 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
## 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
## 99.
## 
## Done.

Quadrat Counting

quadrat <- quadratcount(starbucks_ppp, nx = 10, ny = 10)

plot(
  starbucks_pts,
  cex = 0.5,
  main = "Starbucks Quadrat Counts",
  xlab = "x",
  ylab = "y"
)
plot(quadrat, add = TRUE)

mean(quadrat)
## [1] 10.05333
quadrat.test(quadrat, alternative = "two.sided", method = "MonteCarlo")
## 
##  Conditional Monte Carlo test of CSR using quadrat counts
##  Test statistic: Pearson X2 statistic
## 
## data:  
## X2 = 4821.2, p-value = 0.001
## alternative hypothesis: two.sided
## 
## Quadrats: 75 tiles (irregular windows)
quadrat.test(quadrat, alternative = "regular", method = "MonteCarlo")
## 
##  Conditional Monte Carlo test of CSR using quadrat counts
##  Test statistic: Pearson X2 statistic
## 
## data:  
## X2 = 4821.2, p-value = 1
## alternative hypothesis: regular
## 
## Quadrats: 75 tiles (irregular windows)
quadrat.test(quadrat, alternative = "clustered", method = "MonteCarlo")
## 
##  Conditional Monte Carlo test of CSR using quadrat counts
##  Test statistic: Pearson X2 statistic
## 
## data:  
## X2 = 4821.2, p-value = 5e-04
## alternative hypothesis: clustered
## 
## Quadrats: 75 tiles (irregular windows)

Kernel Density

optim_bw <- bw.diggle(starbucks_ppp, edge = "border")
optim_bw
##    sigma 
## 340.0615
ds1 <- density.ppp(starbucks_ppp, optim_bw)
plot(ds1, main = "Starbucks Density for Optimal Bandwidth")
plot(starbucks_ppp, add = TRUE, cols = "white", pch = 16, size = 0.5)

ds2 <- density.ppp(starbucks_ppp, sigma = 3000)
plot(ds2, main = "Starbucks Density for Bandwidth = 3000")
plot(starbucks_ppp, add = TRUE, cols = "white", pch = 16, size = 0.5)

plot(density(starbucks_ppp), main = "Starbucks Density")
plot(starbucks_ppp, add = TRUE, cols = "white", pch = 16, size = 0.5)

Dirichlet Tesselation

plot(dirichlet(starbucks_ppp), main = "Starbucks Dirichlet Tesselation")

DBSCAN

dbscan1 <- dbscan(starbucks_pts, eps = 500, minPts = 5)
dbscan1
## DBSCAN clustering for 754 objects.
## Parameters: eps = 500, minPts = 5
## Using euclidean distances and borderpoints = TRUE
## The clustering contains 16 cluster(s) and 539 noise points.
## 
##   0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16 
## 539   9   8   5  45   8  16  44  15   6   6   2   5  26   4   5  11 
## 
## Available fields: cluster, eps, minPts, metric, borderPoints
hullplot(starbucks_pts, dbscan1$cluster, main = "DBSCAN Clusters")
## Warning in hullplot(starbucks_pts, dbscan1$cluster, main = "DBSCAN Clusters"):
## Not enough colors. Some colors will be reused.

dbscan2 <- dbscan(starbucks_pts, eps = 2000, minPts = 5)
dbscan2
## DBSCAN clustering for 754 objects.
## Parameters: eps = 2000, minPts = 5
## Using euclidean distances and borderpoints = TRUE
## The clustering contains 9 cluster(s) and 134 noise points.
## 
##   0   1   2   3   4   5   6   7   8   9 
## 134   5 536  20   6  10  15   8   5  15 
## 
## Available fields: cluster, eps, minPts, metric, borderPoints
hullplot(starbucks_pts, dbscan2$cluster, main = "DBSCAN Clusters")
## Warning in hullplot(starbucks_pts, dbscan2$cluster, main = "DBSCAN Clusters"):
## Not enough colors. Some colors will be reused.

dbscan3 <- dbscan(starbucks_pts, eps = 5000, minPts = 5)
dbscan3
## DBSCAN clustering for 754 objects.
## Parameters: eps = 5000, minPts = 5
## Using euclidean distances and borderpoints = TRUE
## The clustering contains 7 cluster(s) and 23 noise points.
## 
##   0   1   2   3   4   5   6   7 
##  23 674   6  10  15   5   6  15 
## 
## Available fields: cluster, eps, minPts, metric, borderPoints
hullplot(starbucks_pts, dbscan3$cluster, main = "DBSCAN Clusters")

dbscan4 <- hdbscan(starbucks_pts, minPts = 5)
dbscan4
## HDBSCAN clustering for 754 objects.
## Parameters: minPts = 5
## The clustering contains 43 cluster(s) and 320 noise points.
## 
##   0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19 
## 320   5   6   6  10  15  14   6   6  10   5  26   8  12  15  10   7  10   6  18 
##  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39 
##   7  10   6   6   5   7  19   8  32  10  10   8   5   7  10  20  15   5  10   5 
##  40  41  42  43 
##  13   7   6   8 
## 
## Available fields: cluster, minPts, coredist, cluster_scores,
##                   membership_prob, outlier_scores, hc
hullplot(starbucks_pts, dbscan4$cluster, main = "HDBSCAN Clusters")
## Warning in hullplot(starbucks_pts, dbscan4$cluster, main = "HDBSCAN Clusters"):
## Not enough colors. Some colors will be reused.

plot(dbscan4, show_flat = TRUE)

Inhomogeneous Poisson Process

model1 <- ppm(starbucks_ppp ~ 1, Poisson())
summary(model1)
## Point process model
## Fitted to data: starbucks_ppp
## Fitting method: maximum likelihood
## Model was fitted analytically
## Call:
## ppm.formula(Q = starbucks_ppp ~ 1, interaction = Poisson())
## Edge correction: "border"
##  [border correction distance r = 0 ]
## --------------------------------------------------------------------------------
## Quadrature scheme (Berman-Turner) = data + dummy + weights
## 
## Data pattern:
## Planar point pattern:  754 points
## Average intensity 1.54e-07 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 10 vertices
## enclosing rectangle: [310427.1, 397312.9] x [3401137, 3499530] units
##                      (86890 x 98390 units)
## Window area = 4894750000 square units
## Fraction of frame area: 0.573
## 
## Dummy quadrature points:
##      64 x 64 grid of dummy points, plus 4 corner points
##      dummy spacing: 1357.589 x 1537.401 units
## 
## Original dummy parameters: =
## Planar point pattern:  2410 points
## Average intensity 4.92e-07 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 10 vertices
## enclosing rectangle: [310427.1, 397312.9] x [3401137, 3499530] units
##                      (86890 x 98390 units)
## Window area = 4894750000 square units
## Fraction of frame area: 0.573
## Quadrature weights:
##      (counting weights based on 64 x 64 array of rectangular tiles)
## All weights:
##  range: [110000, 2090000]    total: 4.89e+09
## Weights on data points:
##  range: [110000, 1040000]    total: 3.95e+08
## Weights on dummy points:
##  range: [110000, 2090000]    total: 4.49e+09
## --------------------------------------------------------------------------------
## FITTED :
## 
## Stationary Poisson process
## 
## ---- Intensity: ----
## 
## 
## Uniform intensity:
## [1] 1.540426e-07
## 
##              Estimate       S.E.   CI95.lo   CI95.hi Ztest      Zval
## log(lambda) -15.68604 0.03641785 -15.75741 -15.61466   *** -430.7238
## 
## ----------- gory details -----
## 
## Fitted regular parameters (theta):
## log(lambda) 
##   -15.68604 
## 
## Fitted exp(theta):
##  log(lambda) 
## 1.540426e-07
AIC(model1)
## [1] 25164.54
model2 <- ppm(starbucks_ppp, ~ log(x) + log(y), Poisson())
summary(model2)
## Point process model
## Fitted to data: starbucks_ppp
## Fitting method: maximum likelihood (Berman-Turner approximation)
## Model was fitted using glm()
## Algorithm converged
## Call:
## ppm.ppp(Q = starbucks_ppp, trend = ~log(x) + log(y), interaction = Poisson())
## Edge correction: "border"
##  [border correction distance r = 0 ]
## --------------------------------------------------------------------------------
## Quadrature scheme (Berman-Turner) = data + dummy + weights
## 
## Data pattern:
## Planar point pattern:  754 points
## Average intensity 1.54e-07 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 10 vertices
## enclosing rectangle: [310427.1, 397312.9] x [3401137, 3499530] units
##                      (86890 x 98390 units)
## Window area = 4894750000 square units
## Fraction of frame area: 0.573
## 
## Dummy quadrature points:
##      64 x 64 grid of dummy points, plus 4 corner points
##      dummy spacing: 1357.589 x 1537.401 units
## 
## Original dummy parameters: =
## Planar point pattern:  2410 points
## Average intensity 4.92e-07 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 10 vertices
## enclosing rectangle: [310427.1, 397312.9] x [3401137, 3499530] units
##                      (86890 x 98390 units)
## Window area = 4894750000 square units
## Fraction of frame area: 0.573
## Quadrature weights:
##      (counting weights based on 64 x 64 array of rectangular tiles)
## All weights:
##  range: [110000, 2090000]    total: 4.89e+09
## Weights on data points:
##  range: [110000, 1040000]    total: 3.95e+08
## Weights on dummy points:
##  range: [110000, 2090000]    total: 4.49e+09
## --------------------------------------------------------------------------------
## FITTED :
## 
## Nonstationary Poisson process
## 
## ---- Intensity: ----
## 
## Log intensity: ~log(x) + log(y)
## 
## Fitted trend coefficients:
##  (Intercept)       log(x)       log(y) 
## -1400.234361     3.253521    89.217040 
## 
##                 Estimate       S.E.      CI95.lo      CI95.hi Ztest       Zval
## (Intercept) -1400.234361 87.7026941 -1572.128483 -1228.340239   *** -15.965694
## log(x)          3.253521  0.7566339     1.770546     4.736497   ***   4.299995
## log(y)         89.217040  5.6957551    78.053565   100.380515   ***  15.663777
## 
## ----------- gory details -----
## 
## Fitted regular parameters (theta):
##  (Intercept)       log(x)       log(y) 
## -1400.234361     3.253521    89.217040 
## 
## Fitted exp(theta):
##  (Intercept)       log(x)       log(y) 
## 0.000000e+00 2.588132e+01 5.577866e+38
AIC(model2)
## [1] 24909.49
plot(model2)

simulation1 <- simulate(model1, nsim = 3)
## Generating 3 simulated patterns ...1, 2, 
## 3.
simulation1[[4]] <- starbucks_ppp
plot(simulation1, main = "")

simulation2 <- simulate(model2, nsim = 3)
## Generating 3 simulated patterns ...1, 2, 
## 3.
simulation2[[4]] <- starbucks_ppp
plot(simulation2, main = "")

Cluster Process

model3 <- kppm(starbucks_ppp, clusters = "Thomas")
summary(model3)
## Stationary cluster point process model
## Fitted to point pattern dataset 'starbucks_ppp'
## Fitted by minimum contrast
##  Summary statistic: K-function
## Minimum contrast fit (object of class "minconfit")
## Model: Thomas process
## Fitted by matching theoretical K function to starbucks_ppp
## 
## Internal parameters fitted by minimum contrast ($par):
##        kappa       sigma2 
## 4.527876e-10 2.078070e+07 
## 
## Fitted cluster parameters:
##        kappa        scale 
## 4.527876e-10 4.558585e+03 
## Mean cluster size:  340.6065 points
## 
## Converged successfully after 127 function evaluations
## 
## Starting values of parameters:
##        kappa       sigma2 
## 1.540426e-07 1.708917e+06 
## Domain of integration: [ 0 , 21720 ]
## Exponents: p= 2, q= 0.25
## 
## ----------- TREND  -----
## Point process model
## Fitted to data: X
## Fitting method: maximum likelihood (Berman-Turner approximation)
## Model was fitted using glm()
## Algorithm converged
## Call:
## ppm.ppp(Q = X, trend = trend, rename.intercept = FALSE, covariates = covariates, 
##     covfunargs = covfunargs, use.gam = use.gam, forcefit = TRUE, 
##     improve.type = ppm.improve.type, improve.args = ppm.improve.args, 
##     nd = nd, eps = eps)
## Edge correction: "border"
##  [border correction distance r = 0 ]
## --------------------------------------------------------------------------------
## Quadrature scheme (Berman-Turner) = data + dummy + weights
## 
## Data pattern:
## Planar point pattern:  754 points
## Average intensity 1.54e-07 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 10 vertices
## enclosing rectangle: [310427.1, 397312.9] x [3401137, 3499530] units
##                      (86890 x 98390 units)
## Window area = 4894750000 square units
## Fraction of frame area: 0.573
## 
## Dummy quadrature points:
##      64 x 64 grid of dummy points, plus 4 corner points
##      dummy spacing: 1357.589 x 1537.401 units
## 
## Original dummy parameters: =
## Planar point pattern:  2410 points
## Average intensity 4.92e-07 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 10 vertices
## enclosing rectangle: [310427.1, 397312.9] x [3401137, 3499530] units
##                      (86890 x 98390 units)
## Window area = 4894750000 square units
## Fraction of frame area: 0.573
## Quadrature weights:
##      (counting weights based on 64 x 64 array of rectangular tiles)
## All weights:
##  range: [110000, 2090000]    total: 4.89e+09
## Weights on data points:
##  range: [110000, 1040000]    total: 3.95e+08
## Weights on dummy points:
##  range: [110000, 2090000]    total: 4.49e+09
## --------------------------------------------------------------------------------
## FITTED :
## 
## Stationary Poisson process
## 
## ---- Intensity: ----
## 
## 
## Uniform intensity:
## [1] 1.542224e-07
## 
##              Estimate       S.E.   CI95.lo   CI95.hi Ztest      Zval
## (Intercept) -15.68487 0.03641785 -15.75625 -15.61349   *** -430.6918
## 
## ----------- gory details -----
## 
## Fitted regular parameters (theta):
## (Intercept) 
##   -15.68487 
## 
## Fitted exp(theta):
##  (Intercept) 
## 1.542224e-07 
## 
## ----------- CLUSTER  -----------
## Model: Thomas process
## 
## Fitted cluster parameters:
##        kappa        scale 
## 4.527876e-10 4.558585e+03 
## Mean cluster size:  340.6065 points
## 
## Final standard error and CI
## (allowing for correlation of cluster process):
##              Estimate      S.E.   CI95.lo   CI95.hi Ztest      Zval
## (Intercept) -15.68487 0.6235393 -16.90698 -14.46276   *** -25.15458
## 
## ----------- cluster strength indices ----------
## Sibling probability 0.8942624
## Count overdispersion index (on original window): 298.773
## Cluster strength: 8.457373
## 
## Spatial persistence index (over window): 0
## 
## Bound on distance from Poisson process (over window): 1
##   = min (1, 1509.76, 5076500, 4819384, 2195.309)
## 
## Bound on distance from MIXED Poisson process (over window): 1
## 
## Intensity of parents of nonempty clusters: 4.527876e-10
## Mean number of offspring in a nonempty cluster: 340.6065
## Intensity of parents of clusters of more than one offspring point: 4.527876e-10
## Ratio of parents to parents-plus-offspring: 0.002927345 (where 1 = Poisson 
## process)
## Probability that a typical point belongs to a nontrivial cluster: 1
model4 <- kppm(starbucks_ppp, clusters = "LGCP")
summary(model4)
## Stationary Cox point process model
## Fitted to point pattern dataset 'starbucks_ppp'
## Fitted by minimum contrast
##  Summary statistic: K-function
## Minimum contrast fit (object of class "minconfit")
## Model: Log-Gaussian Cox process
##   Covariance model: exponential
## Fitted by matching theoretical K function to starbucks_ppp
## 
## Internal parameters fitted by minimum contrast ($par):
##       sigma2        alpha 
##     2.806734 11590.392857 
## 
## Fitted covariance parameters:
##          var        scale 
##     2.806734 11590.392857 
## Fitted mean of log of random intensity:  -17.08824
## 
## Converged successfully after 59 function evaluations
## 
## Starting values of parameters:
##   sigma2    alpha 
##    1.000 1307.255 
## Domain of integration: [ 0 , 21720 ]
## Exponents: p= 2, q= 0.25
## 
## ----------- TREND  -----
## Point process model
## Fitted to data: X
## Fitting method: maximum likelihood (Berman-Turner approximation)
## Model was fitted using glm()
## Algorithm converged
## Call:
## ppm.ppp(Q = X, trend = trend, rename.intercept = FALSE, covariates = covariates, 
##     covfunargs = covfunargs, use.gam = use.gam, forcefit = TRUE, 
##     improve.type = ppm.improve.type, improve.args = ppm.improve.args, 
##     nd = nd, eps = eps)
## Edge correction: "border"
##  [border correction distance r = 0 ]
## --------------------------------------------------------------------------------
## Quadrature scheme (Berman-Turner) = data + dummy + weights
## 
## Data pattern:
## Planar point pattern:  754 points
## Average intensity 1.54e-07 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 10 vertices
## enclosing rectangle: [310427.1, 397312.9] x [3401137, 3499530] units
##                      (86890 x 98390 units)
## Window area = 4894750000 square units
## Fraction of frame area: 0.573
## 
## Dummy quadrature points:
##      64 x 64 grid of dummy points, plus 4 corner points
##      dummy spacing: 1357.589 x 1537.401 units
## 
## Original dummy parameters: =
## Planar point pattern:  2410 points
## Average intensity 4.92e-07 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 10 vertices
## enclosing rectangle: [310427.1, 397312.9] x [3401137, 3499530] units
##                      (86890 x 98390 units)
## Window area = 4894750000 square units
## Fraction of frame area: 0.573
## Quadrature weights:
##      (counting weights based on 64 x 64 array of rectangular tiles)
## All weights:
##  range: [110000, 2090000]    total: 4.89e+09
## Weights on data points:
##  range: [110000, 1040000]    total: 3.95e+08
## Weights on dummy points:
##  range: [110000, 2090000]    total: 4.49e+09
## --------------------------------------------------------------------------------
## FITTED :
## 
## Stationary Poisson process
## 
## ---- Intensity: ----
## 
## 
## Uniform intensity:
## [1] 1.542224e-07
## 
##              Estimate       S.E.   CI95.lo   CI95.hi Ztest      Zval
## (Intercept) -15.68487 0.03641785 -15.75625 -15.61349   *** -430.6918
## 
## ----------- gory details -----
## 
## Fitted regular parameters (theta):
## (Intercept) 
##   -15.68487 
## 
## Fitted exp(theta):
##  (Intercept) 
## 1.542224e-07 
## 
## ----------- COX  -----------
## Model: log-Gaussian Cox process
## 
##  Covariance model: exponential
## Fitted covariance parameters:
##          var        scale 
##     2.806734 11590.392857 
## Fitted mean of log of random intensity: -17.08824
## 
## Final standard error and CI
## (allowing for correlation of Cox process):
##              Estimate      S.E.  CI95.lo   CI95.hi Ztest      Zval
## (Intercept) -15.68487 0.7361998 -17.1278 -14.24195   *** -21.30518
simulation3 <- simulate(model3, nsim = 3)
simulation3[[4]] <- starbucks_ppp
plot(simulation3, main = "")

simulation4 <- simulate(model4, nsim = 3)
simulation4[[4]] <- starbucks_ppp
plot(simulation4, main = "")